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Abstract. Dynamic NURBS, also called D-NURBS, is a known dynamic version of the 
nonuniform rational B -spline (NURBS) which integrates free-form shape representation 
and a physically-based model in a unified framework. More recently, computer aided de- 
sign (CAD) and finite element (FEM) community realized the need to unify CAD and 
FEM descriptions which motivates a review of D-NURBS concepts. Therefore, in this 
paper we describe D-NURBS theory in the context of ID shape deformations. We start 
with a revision of NURBS for parametric representation of curve spaces. Then, the La- 
grangian mechanics is introduced in order to complete the theoretical background. Next, 
the D-NURBS framework for ID curve spaces is presented as well as some details about 
constraints and numerical implementations. In the experimental results, we focus on pa- 
, rameters choice and computational cost. 
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\0 ■ 1 Introduction 

2 In the context of animation of soft objects every engine is composed by three linked parts: the geometric model, 

^ dynamic model and rendering module. The former can be realized in the context of parametric frameworks like 

^—^ . nonuniform rational B -spline (NURBS) ||Piegl and Tiller 1997[ IFarin 199711. The dynamic model needs physic 



X 



models that incorporate dynamic quantities like velocity, mass and force distributions, into an evolution equation 
that governs the shape deformation HErleben et al. 20051 . The latter includes global/local illumination techniques 



to generate the scene with the desired realism [?]. In this work we focus only on the first two components. 



Non-uniform Rational B-spline (NURBS) is a mathematical framework commonly used for generating and 
representing curves, surfaces and volumes | Piegl and Tiller 1997) . It offers an unified mathematical basis to 



describe analytic and free-form shapes with great flexibility and precision. NURBS became a standard for CAD 
(Computer Aided Design) systems due to its excellent mathematical, numeric and algorithmic properties. NURBS 
are built from the B-spline function basis and a NURBS curve is a composition of NURBS functions, a set of 
control points {pi, P2, • • ■, Pn} C and a weight vector {wi,W2, ■ ■ -^Wn)- The control points and the weights 
compose the degrees of freedom of the NURBS curve. 

For computer graphics applications, the dynamic model in general is based on classical mechanics which 
is concerned with physical laws to describe the behavior of a macroscopic system under the action of forces 
IIDeusen et al. 20041 . For instance, when considering a particle in the 3D space under the action of gravity, we 
can take its position vector along the time t, which in cartesian coordinates is given by {x{t),y{t),z{t)), and 
use the Newton's laws to get the governing equation written in terms of the cartesian coordinates and the time 
t. In a more general situation, the instantaneous configuration of a system may be described by the values of n 
generalized coordinates {pi,P2, ■ • ■,Pn)- So, we need a methodology to write the evolution equation of the system 
in terms of the generalized coordinates. 

The Lagrangian formulation of mechanics is a framework to address this issue HGoldstein 19811 . It is a 
variational formulation of mechanics based on the integral Hamilton's Principle which states that the motion of 
the system between times ti and t2 is derivable from the solution of a variational problem HGoldstein 198 1 II . 



The corresponding Lagrange's equations allow to write the evolution of the system in term of the generalized 
coordinates. That is what we need to link the geometric model of NURBS and the dynamic model: we can use 
the control points and the weights as generalized coordinates to describe the physical system. Therefore, we get 
an approach that integrates shape representation and a dynamic model in a unified framework called D-NURBS 



in the literature |Terzopoulos and Qin 1994 Qin and Terzopoulos 1996 1. 



Continuous systems, like an elastic curve, have infinite degrees of freedom which difficult its description for 
both the geometric and dynamic aspects. In mathematical terms, we are dealing with infinite basis functions, 
may be uncountable. One possibility to simplify the problem is to consider finite dimensional representation 
with enough flexibility in order to represent the solution with the desired precision. In the context of mechanical 
systems the Finite Element Method (FEM) is the traditional way to perform this task. However, as pointed out 
in liCottrell et al. 2009 J . NURBS framework can be also considered. That is way geometric modeling and FEM 
community realized the need to unify CAD and FEM descriptions which motivates our review of NURBS and 
D-NURBS concepts. 

So, in section |2] we start with an objective review of B-splines functions in order to set up the background 
for NURBS development. Next, in section |3] we describe the Lagrangian mechanics framework in the presence 
of constraints and a generalized potential for dissipation forces. Then, section |4] considers the D-NURBS model 



following the presentation given in [Terzopoulos and Qin 1994 1. We present details of the evolution equation 



generation, constraints introduction and numerical aspects. For simplicity, we focus on curve spaces but the 
theory can be straightforward generalize for surfaces and volumes. In the experimental results (section |5]), we 
consider a set up for a linear mass distribution with fixed end points. We discuss the influence of parameters 
choice, effects of NURBS weights and computational cost. The conclusions and further works are presented in 
section [6l The appendices |A] and |B] give some details about specific terms of the D-NURBS governing equation. 

2 NURBS: Nonuniform Rational B-spline 

The spline framework is the starting point for NURBS development. A polynomial spline of order k (degree k — 1) 
is a piecewise polynomial function of order k with continuity of derivative of order k — 2 at the common joints 
between segments, which are called patches [Rogers and Adams 1976 IPersiano 19961 . 



Therefore, the spline space is a functional space composed by piecewise polynomial functions with the prop- 
erty already stated. A fundamental element in the spline theory is the knot vector which defines the end points 
of the patches of the spline function. Given the order k and a knot vector v = (uo,ui, ■ ■ -^Un), we can denote 
the space of polynomial splines of order k with domain in the range [uq, n„] as S^{uq,ui, ■ ■ -jUn)- The Figured] 
shows some elements of this set when k = 1, 2, 3. 



We can show that the S'^^uq, ui, • • •, u^) is a vector space of dimension n — k + 1 UPersiano 19961 . The main 
point in the spline theory is to construct a basis for this space. From the functional analysis viewpoint the space 
properties are invariant respect to the basis choice. However, for computer graphics aspects it is important that 
every tool and algorithm generated has an intuitive geometric and visual interpretation with local control of the 
target objects. The B-spline basis attend these requirements. 

Following traditional texts in this area UFarin 19971 [Rogers and Adams 1976[ we perform a recursive defini- 



tion of the B-spline basis. So, let us consider the S^{uq, ui, n^); that means, the space of piecewise polynomial 
functions of order k = I (degree 0) which are just piecewise constant functions, like the one presented on Figure 
□for k = l. 

A basis for this space is in fact the first B-spline basis in our recursive scheme, which is defined as follows: 

„ .X fl, ifui<u< Ui+i 
I U, otherwise 

for i = 0, 1, • • •, n — 1. 

Now, let us consider the space S'^{uq,ui, ■ ■ ■,Un)] that is, the space of piecewise polynomial functions of 
order k = 2 (degree k — 1 = 1) which are just piecewise linear functions with continuity of derivative of order 




k=2 



k=1 



Figure 1: Polynomial spline examples with knot vector v = (^0)^11*21^31*4, is)- 



k — 2 = (see Figure [U. We are supposing that mq < "Ui < • • • < Un- We already know that this is a vector 
space with dimension n — k + l = n — 1. Besides, when integrating polynomial functions of order k we get 
again polynomial functions but with order k + I. Also, we want that the support of the functions i?j 2 would be 
as small as possible (for local geometric control) and that they have continuity of derivative of order — 2 = at 
the common joints between patches. The following functions fulfill these requirements: 

Bi,2{u)= ■ ■ ] ds, 1 = 0,1,- ■ ■,n-2. (2) 

J_oo\Ui+l-Ui Ui+2-Ui+lJ 

By repeating the above arguments, we can show that the following recursive scheme will generate a basis 
= {Bi^k, z = 0, 1, • • ■, n — /c — 1} for splines / : [uq, Un] —> 5ft such that / (uq) = / (un) = 0: 

„ / N (n - Uj) Bj^k-i ju) [uj+k - u) Bj+i^k-i ju) . „ , , 

Bi,k{U) = \ , I = 0,1,- ■ -,71 - k (3) 

Ui+k^l — Ui Ui+k — Ui+l 

where k = 2,3, ■ ■ ■ and the i?j 1 (n) is given by expression ([T}. The Figure |2]pictures the obtained basis for k = 2. 



So, in the above development, the span of B^ is in fact a subspace of S^{uq,ui, ■ ■ ■,Un) once B'' can only 
generate functions with support in the interval {uQ,Un) as we already observed above. However, we can cover all 
the spline space by considering more general knot vectors. In fact, the knot vector has a significant influence in 
the spline basis generated. In general, it is used three types of knot vectors: uniform, open uniform (or just open) 
and nonuniform. 

Uniform knot vectors satisfies Uj+i — Ui = Au = const., for i = 1,2, ...,n. Uniform knot vectors yield 
periodic uniform basis functions, like the one presented in Figured that means: 

Bi^k {u) = A-i,fe {u - Au) = Bi^i^k {u + An) . 



An open uniform knot vector has also the property ui+i — Ui = Au for internal knots but it has multiplicity 
of knot values at the ends equal to the order k of the B-spline functions. For instance: 




Figure 2: B-spline of order k = 2 with knot vector v = {uo,ui,U2, us,U4, u^). 



v = (0,0,1,2,3,4,5,5), if k = 2, 



V = (0,0,0,0.1,0.2,0.3,0.4,0.5,0.5,0.5), if k = 3. 

These kind of knot vectors may yield more general B-spline basis that can generate functions that are not 
null at the ends of the knot vector, as we can visualize in Figure |4] 



Finally, nonuniform knot vectors may have either unequally spaced (ui+i — Ui = Aui) and/or multiple knot 
values at the ends or even for the internal knots. 

The B -splines generated by open (uniform or nonuniform) knot vectors have important properties [ Piegl and Tiller 1997] . 

1- Bi^k (u) > Vu. 

2. Sj fc (u) = if u is outside the interval G [uj, tij+fc+i). 

3. Partition of unity: Y^'I^q Bi^k (u) = 1. 

Once defined the basis for the spline space, we can consider curve spaces in generated through B-splines. So, 
let us take a set of pi, i = 0,1,2, ■ ■ ■,n — k points in 5)?^ and the vector-valued function given by: 

n—k 

c{u) = Y^ PiB^^k (u) . (4) 

i=0 

This function defines a curve of class C'^~^ in 9?^, which is called a spline curve. The points Pi are called 
control points and the corresponding polygon is the defining polygon. Important properties about these curves 
are: 

1. End points interpolation: in the case of open knot vector we have c (mq) = Po and c (u„) = p„. 

2. Affine Invariance: If ip (r) = Ar+v is an affine transformation then tp (c (u)) = J^IZ^ Bi^k {u) ip (pj). 




Figure 3: B-splines examples of order k = 2 with uniform knot vector v = (0, 1, 2, 3, 4, 5). 

3. Strong convex hull property: the curve belongs to the convex hull of its control polygon. 
A rational B -spline curve is the projection of a polynomial B -spline curve defined in the four-dimensional homo- 



geneous coordinate space back into the three-dimensional physical space [Rogers and Adams 1976J . Therefore, 



if we represent the control points in the four-dimensional homogeneous coordinate space we obtain: 

Pi=i )> i = 0, 1,2,- • -,71 - A;, 



and applying expression © we get a spline curve in the four-dimensional homogeneous space: 

n—k 



c\u) = 

1=0 

By projection in the three-dimensional space we obtain the rational curve 



En-k D / \ n—k 

c[n) = —^^ — — — = 2^PiiV,,fc(tx), (6) 



^n—k D / \ / -I 

where Ni^^ the rational B-spline functions given by: 



AT ( \ WjBj^k {U) 

J^ik{u) = ; — • (7) 

If the B-splines in expression |3]are generated by nonuniform knot vectors then the functions i?j ^ in expression 
([7]) are named nonuniform rational B-splines (NURBS) and the curve defined by expression (|6j is a NURBS curve. 

B-splines can be enriched without modifying the underlying geometry and parameterization through the 
mechanisms that are called refinements. The most common mechanisms are knot insertion and degree elevation 
IIFarin 19971 |Piegl and Tiller 1997| . 



3 Lagrangian Mechanics 

Let us consider a physical system whose instantaneous configuration may be described by the values of n general- 
ized coordinates {pi,P2,- ■ -^Pn) which can be considered as a point in a n — dimensional Cartesian hyperspace 




Figure 4: B-splines functions of order k = 2 with open uniform knot vector v = (0, 0, 1, 2, 3, 4, 5, 5). 

known as configuration space. As time goes on from a time ti to a time t2, the system changes its configuration 
due to internal and external forces. Therefore, the evolution of the system can be seem as a continuous path, or 
curve, p(t) in the, configuration space, parameterized through the time t. 

The Hamilton's Principle gives a methodology to write the evolution equation of the system in terms of the 

generalized coordinates and time t. It states that if for a mechanical systems with kinetic energy T = T ^p^ , 

where p = dp/dt, all force fields are derivable from a scalar potential V = V {p, t) then the motion of the system 
from time ti to time t2 is such that the line integral: 

I = L(p,p,tjdt, (8) 

where L (^p,p,t^ = T ^p^ — V (p, t) , has a stationary value for the correct path of the motion HGoldstein 19811 . 

The function L is named the Lagrangian of the system and we can apply traditional techniques of the varia- 
tional calculus to show that the correct path must satisfies: 

d ( dL\ dL 

or, in a compact form: 

'dL\ dL 

— - — = 0, i = l,2,---,n, (10) 
dp J 9p 

which are the Lagrange equations of motion HGoldstein 19811 . 

We can introduce dissipation forces in the Hamilton' principle by adding a velocity-dependent term in the 
scalar potential of the system. So, let us consider the general form for the Lagrangian: 



L p,p,i =r p,p - [/(p) + F p,p , (11) 



where, like before, T is the kinetic energy but now possibly dependent from both p and p, [/ is the potential 
related to the conservative forces and F is a velocity-dependent potential to account for dissipative effects. 



So, substituting this expression in the Euler-Lagrange equations (fTOl) renders: 



d dT dF\ /dT dU dF 

dt \dp dp J V^P ^P 

In general, mechanical systems undergoes effects of internal and external forces. Therefore, it is useful to 
decompose the potential U (p) into two terms named Eint e E^xt^ which will account for the internal and external 
forces, respectively: 

U (p) = E,nt (p) + E,xt (p) . (13) 
By substituting expression ( fT3] ) into the equations (fT2l ). we get: 

d fdT dF\ ^ Z dEinA _ dEext ^ fdT dF\ ^^^^ 
dt\dp dp) \ 9p J dp \dp dp J' 

which gives the general form of Euler-Lagrange equations. 




3.1 Lagrange Equations with Constraints 

Now, let us extend the Hamilton' principle in order to cover constraints. We focus on holonomic constraints; or 
holonomic system, for which the constraints may be expressed by: 



fl{Pl,P2,- ■ ■,Pn,t) = 0, 
f2{pi,P2,- ■ ■,Pn,t) = 0, 

(15) 

fm{pi,P2,- ■ ■,Pn,t) = 0, 

where //, / = 1, 2, • • •, m, is a general expression connecting the generalized coordinates. In this case, we can 
take the differential dfi: 

dfi = y2^dpk + ^dt = 0, / = l,2,---,m. (16) 

^ 9Pk ot 

If we consider dt = and replace dpk by the corresponding virtual displacement Spk we can rewrite expres- 
sion ( fT6l ) as: 

n 

^aikSpk = 0, l = l,2,---,m, (17) 

fc=i 

where aik = dfi/dpk - Expression ([17] ) implies a dependence between the virtual displacements 5pk- In order to 
reduce the number of virtual displacements to only independent ones we can use Lagrange multipliers Ai, A2, • • 
•, Am. So, we can put together the equations ([TT] ) using the expression: 

/ J]J^Aia;fc(5pfc = 0. (18) 

k=i 1=1 

Therefore, by assuming that the Hamilton's principle holds for holonomic systems we can incorporate ex- 
pression ( fTSl ) in the variational technique used to get Lagrange equations ^ and to obtain: 



rt2 ( dL d ( dL\ \ 



5pk = 0. (19) 



We shall remember that the virtual displacements 5qk are connected by the m equations (fTTl) . Besides, the 
Lagrange multipliers Ai , A2 , • ■ ■ ) remains at our disposal. So, let us suppose that we can choose these multipliers 
such that: 



dL d ( dL\ ^ ^ ^ , 

— — + > ^ Xiaik = 0, k = n-m + l,- ■ -,171. 



(20) 



By substituting this expression in the integral ([T9l ) we render: 

„^2 n—m ' ' ^ 



(21) 



Once we have m constraint equations in ([TtT i the only virtual displacements 6pk involved in expression (I2TI 1 
are the independent ones. Therefore, it follows that: 



dL d dL\ ^ ^ 
opk az \dpkj 1=1 



,n — m. 



(22) 



Expressions (1201 ) and (122] ) give the complete set of Lagrange's equations for holonomic systems. However, 
the expressions involves n + m unknowns, namely the n coordinates p^ and the m multipliers A;. So, we must 
add to the final result the constraints give by expression ([T6] ). 

Therefore, by putting together expressions (|20) . (I22t and ([T?] ) we find that the desired solution must satisfies 
the equations: 



dL d dL\ ^ ^ 
upi ai \dpi / 



fl{Pl,P2,- ■ ■,Pn,t) =0,, 1 = 1,2, 



, m. 



(23) 
(24) 



4 D-NURBS Formulation 

The idea is to submit an initial NURBS curve, given by expression (O to a Newtonian dynamics generated by an 
external potential, internal (elastic) and dissipation forces. Therefore, a natural way to parameterize the evolution 
of the curve along the time is: 



c{u,t) 



(25) 



So, the control points pi{t) and the weights Wi{t) becomes time-dependent while the rational functions 
(O remain u-dependent only. Therefore, the control points and the weights become the degrees of freedom 
of the system evolution; and so, they compose the generalized coordinates which are concatenated as follows 



I Terzopoulos and Qin 1994) : 

p{t)=[{plwo) {pIw,) ... {plwn)f 
where we have the control points vector and the weights vector specified, respectively, by: 

p^{t) = [wo wi ... «;„f€K("+i). 
A fundamental element in the D-NURBS development is the associated Jacobian, defined as follows: 



J 



Bo {u, p) 



9c 



Bniu,p),^ 



G 



Sj^3x4(n+l) 



(26) 

(27) 
(28) 

(29) 



where: 



Bi {u, p) = 
with (see expression |7]): 



i9c <9c 8c 

dpix dpiy dpi; 



Ni^k {u, Pii 










and: 



9c YJj=o{Pi{t) - Vj{t))wjBi^k {u) Bj^k {u) 



dwi 



(30) 



(31) 



(32) 



We shall observe that Bi {u, p) G 3?^^^ and ^ G 3?'^, com i = 0,. . . ,n and consequently J G Sfj3x4(n+i) 



We can concatenate the Bi's and according to the following matrices: 

B = [Bo{u, p), Bi{u, p), • • •, Bn{u, p)] G K3x3(n+1) 
5c 5c Oc 



ly 



G 



Sj^3x(n+l) 



(33) 
(34) 



c?Tz;o dwi dwn 

The advantages of defining the matrices J, W and the vectors pf, and pai becomes clear by observing that: 



Jp 



POx 
POy 
POz 
Wo 
Plx 
Ply 
Plz 
Wl 



YA=oWi{t)Bi^k{u)Vi{t) ^ 



Eto(P*W - Pj(.t))wj{t)B,^k{u)Bj,k{^) 



\ 



But, with a simple algebra we can show that: 



Bpb + Wp 

w 

Wp^ = 0. 



Pnx 
Pny 
Pnz 



Wi{t) 



(35) 



(36) 



Therefore, 



Jp = Spb. 

However, by remembering expression (l25T l it is clear that: 



(37) 



c{u,t) = Bpb. 
Henceforth, from expressions (l37l l and (l38T i we get that: 



(38) 



c(u,p) = Jp. 
Other important properties that can be easily proved are: 

dJ , , 

^ ■•>(" = »■ 

dc{u, p) J dp 
di ~ ' lit' 



(39) 

(40) 
(41) 



The next step is to compute the kinetic and (generalized) potential terms to be inserted in the Lagrangian 
given by expression ([TTI i. 

4.1 Kinetic Energy T 

In this work we focus on the D-NURBS formulation for a continuous parametric curve subject to a force field. So, 
we shall consider a (constant) linear mass density distribution fi. Therefore, the kinetic energy is computed by: 



T 



1 





dc 


1 u 


lit 



du, 



(42) 



where ^ is the curve velocity. By applying expression (1411 we observe that: 



dc 

lit 



{Jpf ■ (Jp) 



So, if we insert expression (l43T i into kinetic energy (|42] | we obtain: 

T=^l^fx{Jpf .{Jp)du. 



which becomes: 



T = ^ / np^J^Jpdu, 



Once p does not depend on the parameter u, we can rewrite expression 

T = ^p^Mp, 



as: 



(43) 



(44) 



(45) 



(46) 



M = M(p) = / flJ^Jdu G Sft4(n+I)x4(™+1)^ 
J u 



(47) 



where: 

is called the mass matrix. 
4.2 Energy Dissipation F 

Formally, the idea is to consider a velocity-dependent potential F such that, when introduced in the Euler- 
Lagrange equations ([141 ) generates a velocity-dependent dissipative force. In order to perform this task let us 
suppose that F satisfies: 



dF _ 




dc 






lii 



du 



F{t) 



t=0 Ju 



7 



dc 
lit 



dudt, 



(48) 



where the constant 7 is the the damping density. By performing an analogous development of section 14.11 we 
obtain: 



dF 



1 



where D € S)[^4(n+i)x4(n+i)^ damping matrix, is computed by: 



D = D{p) = I -fJ^Jdu. 

J u 



(49) 



(50) 



4.3 Potential for Conservative Forces 

The internal and external conservative forces are introduced in the D-NURBS Lagrangian through the potentials 



Eint and Eext, respectively. We compute the former by using the thin-plate model | Terzopoulos and Fleischer 1988 | 



Eint (p) 



a 



dc 



du 



+ /3 



d^c 



dv? 



du, 



(51) 



where a is the elasticity and /3 the rigidity parameter of the curve. Using the expression (|39l ) and the fact that the 
generalize coordinates vector p does not depends on the parameter u (see expression (|26l )) we can show that: 



(52) 



(53) 



dc d 

du du 

Obviously the same is true for the second derivative respect to the parameter u . Therefore: 

Eint (P) = ^ y (" (-^"P)^ (-^"P) + /5 (JuuP)^ (JuuP)^ du. 

Once (JuP)^ (JuP) = P^Ju^uP it follows: 

Eint (P) = ^ y {oip'^JuJuP + f^p'^ JuuJuuP) du, 

and, consequently: 

Eint (P) = ^P^^P' 

where the matrix K = K{p) € S)fj4(n+i)x4(7i+i)^ named the stiffness matrix, is given by: 

K{p) = / {aJ^Ju + PJuuJuu) du. 

J u 

The external potential E^xt generates the force fields, like gravity, that act on the system. According to 
expression ([14] ). they are computed by the gradient of the potential E(.xt respect to the generalized coordinates: 



(54) 



(55) 



(56) 



dEext 1 / dEext dEext QE^xt QE^xt dEext dEext dEext QE^xt dEext dEext dEext dEext 



dp 2\ dpox ' dpoy ' dpoz ' Owq ' dpix ' dpiy ' dpiz ' dwi ' ' dpnx ' dpny ' dpm ' dwn 



(57) 



4.4 Euler-Lagrange Equations for D-NURBS 



Now, we insert the kinetic energy and potentials just computed in the Euler-Lagrange equations given by expres- 
sion (fT4l ). Besides, we must observe that the matrices M, D and K are all symmetric and for a quadratic form 
g = p^Ap with A symmetric we have ^ = 2Ap. Therefore: 



dt 



d / dT 



dp 



i ( I (^P^^P) ) = i (^2Mp) = i (Mp) = Afp + Mp 



d / 



f = |;(ip^Mp) = i(p)''^P 



9p 9p 

_ 1 r* 



^(ip^Kp) = l2i^p + i 



p^fp 



T 



i^p+i 



P^fP 



By substituting these expressions in the Euler-Lagrange equation: 



d IdT dF\ (dEint 
dt V^p 9py V 5p 



/ar_ ap 

(?p V 9p 9p 



we get: 
Mp + Mp 



Dp +i^p+ 



1 rj.dK 

2P 



ext 



5p 



+ 



1 /■\T dM 



(58) 



dp 



2 



p^— pd-u ) dt, 
t=o \Ju dp ) 

(59) 



dD 



which can be rewritten as follows by just re-arranging the terms: 



Mp + Dp + i^p 



dE, 



ext 



dp 



+ 



r dM 



dp 



P 



■Mp- 



1 

2^ 



+ ■ 



t=0 



p"^— pdn ) dt. (60) 
ap 



However, in the Appendices A and B we shown that: 



lp = Mp--(p) — p, 



(61) 



1 T^dK 
2P 



0, 



where: 



/iJ Jdii. 



(62) 



(63) 



Therefore, if we neglect the effects of the last integral term, we can finally write the governing equation for 
D-NURBS as: 

d F 

Mp + Dp + i^p = - /p, (64) 

dp 



where the matrix / is computed by equation (I63I) . 



4.5 External Forces 

In this section we consider an external potential which in cartesian coordinates has the general form: 



E, 



ext 



P (x, y, z) du 



where P {x, y, z) is a potential density function. 
So, 



Eext{p)= / P{x{p,u) ,y{p,u) ,z{p,u))du 

J u 



(65) 



(66) 



In cartesian coordinates, the external force is given by: 



ext= 



dP 

"57 



du. 



(67) 



However, we must write the external force respect to the generalized coordinates, following the expression 
57l) . For instance, let us consider the term: 



ext 



From the Chain-Rule: 



dpox 

dEext 



d 



■P{x (p,u) ,y (p,m) ,z(p,u)) 



'Ox 



But, from the expression (|25] ) we observe that: 

dy 



Opox 

dP dx dP dy dP dz 
\ ^ -I 

dx dpox dy dpox 
dz 



du 



dz dpox 



du. 



Therefore: 

Analogously we can find: 

dEext 



On the other hand: 
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dw( 



P{x{p,u) ,y{p,u) ,z{p,u)) 



du, 



and so: 



dwn 



dP dx dP dy ^ dP dz 



dx dwo dy dwo dz dwo 



du — 



dx dy dz 
dwo ' dwo ' dwo 



dP dP dP 

dx' dy' dz 



du 



(68) 
(69) 

(70) 
(71) 

(72) 
(73) 

(74) 



So, by using expression ( [30l ). we find that the above results can be grouped in the following matricial expres- 



sion: 
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(75) 



Generalizing for z = 0..n we have: 
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(76) 



where / (x, y, z) is the external force field density defined by the gradient of the potential density P in cartesian 
coordinates {x, y, z) . Therefore, the external force field in the generalized coordinates is given by: 

fp{p)= J'^f{x{p,u),y{p,u),z{p,u))du. (77) 

J u 

Particularly, in the case of the gravitational potential for a particle we have: 

E = -mgy, (78) 

where m is the mass particle, g is the gravitational field intensity and y gives the particle position (its height) in 
the vertical axes. For a ID continuous system, a curve in the two-dimensional Euclidean space, the gravitational 
potential can be computed by a generahzation of expression (178] ) given by: 



Eext = - I fJ-gydu (79) 

where /i is the Unear mass density (constant), like before. Therefore, the potential density is: 

P{x,y) = -^lgy. (80) 

and the force field density is given by: 

VP (x, y) = -^^g{^, = -f,g{0, if (81) 

So, according to expression dTTT l. the external force field /p is: 

/p (?) = -/" f^J'^g (l)du = -l'fiJ^(^]du. (82) 



4.6 D-NURBS with Constraints 



In the case of linear constraints, equations ([TST i become: 



+ d = 0, (83) 

where A G K™^", with m < n, is a constant matrix and d € 9f?™ is a constant vector. In this case, we can choose 
a set of n — m, say q = {qi,q2, ■ ■ -^qn-m) independent variables and explicitly write the m remaining ones as a 
function of the q vector, which will be the new generalized coordinates. In fact, if we write equation (1831 ) in the 
form: 
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or simply: 

Giq + G2q = -d 

where Gi € SR'"^'" and G2 G s)fjmx (n-m) ^j^g ^j.^^ ^j^^ second matrices of expression (l84l and q = 
{Pi,P2, ■ ■ ■,Pm) , q = {pm+i,Pm+2, ' ' ■,Pn) ■ Then, by supposing that pi, ■ ■ ■,pm can be choosen such that 
Gi is non-singular, we have: 

q=-Gr^G2q-GrM. (85) 



Let / G sfi{n-m)xin-m) f^^^ jj^g observation that G^^G2 G K"^x("-'™) and using expression ^ it is 
clear that the matrix: 



G 



-G^'G2 
I 



allows to write: 



(86) 



p = Gq+do 

with do = [G];^d,of G K^^^ 

The equations (1231 ) can be written in compact form as: 



a TP 

Mp + Dp + i^p+— £^ + /p 1 . 



5p 



So, using expressions ([87] ) we can observe that: 



p = Gq, p = Gq. 



(87) 



(88) 



(89) 



Therefore, by substituting expressions (ISTT i and (|89] l in equation (1881 ) and using the fac that A = [Gi G2] 
we obtain: 



Gl 



A 



MGq + Z)Gq + K (Gq+do) +^^|? + /Gq 



(?p 5q 

If we multiply both sides by G^, where the matrix G is defined in expression (l86l ) we obtain: 



(Gr'G 



Gl 



dE, 



X = -G' X ( MGq + DGq + K (Gq+do) +^^^ + ^^q 



= -G' MGq - G^Z)Gq - G' K (Gq+do) - G' 



T dEext 



dp 



G^IGq 



G^MGn + G^DGq + G^KGq_ = -G^— ^ - G^/Gq - G^i-Cdo. 

op 

If we name: 

Mg = G^MG; Dg = G^DG; Kg = G^KG; fg = -G^^^'~^ 



dp ' " 

then, we get the D-NURBS evolution equation subject to the linear constraints given by: 

Mgq + Dgq + K^q = /, - /,q - G^i^do. 



; In = G'IG, 



(90) 



4.7 Numerical Implementation 

The equation ( [64l ). as well as its constrained counterpart in expresson (|90l ). does not have in general analytical 
solution and so we have to use a numerical approach to solve it with the desired precision. The equation (l64l ) is a 
second order ordinary differential equation. Besides, it is important to observe that the matrices M, D, K depends 
on the integration of products of the rational B -spline functions dT) and their derivatives of first and second order 
respect to the variable u. 

Therefore, the numerical solution of expression (l64T i can be performed by finite difference methods (FDM) 
in time. Besides, we need a numerical scheme for computing the integrals, as described next. 



4.7.1 Matrices Computation 



The matrices M, D, K that appears in D-NURBS evolution equation are given by expressions (I47l).(l50l). and (l56l) . 
respectively. They involve derivatives of zero, first and second order of J respect to the variable u. For instance, 
for matrix K = (%) G K4(n+i)x4(n+i) behave: 



where: 




with j j means the coUum i of the Jacobian J. 

The computation of each term in the summation in expression (|9T| ) can be performed by Gauss quadrature 
I Chapra and Canale 2009) . An analogous scheme can be used to compute the other matrices. 

In our implementation we have developed a numerical approach based on isogeometric analysis follow- 
ing the recipe of HCottrell et al. 2009]| . Our implementation avoids the cost of assembling the global matrices 
M, D and K. For this, we calculate the matrices of each element individually, where the elements are con- 
structed by partitioning the knots vector. Figure |5] shows building elements (el, e2, e3, e4) from open knots 
vector V = (0, 0, 0, 0.25, 0.50, 0.75, 1, 1, 1). Here we will have six control points (n — fc + l = 8 — 3 + 1 = 6) 
which according MCottrell et al. 20091 will be distributed over the elements by following expression 

E = {el = {pi,p2,P3) ,e2 = (p2,P3,P4) , e3 = {p3,P4,P5) = {p4,P5,P6)} (93) 



v=(0,0,0,0.25,0.50,0.75,l,l,l) 




Figure 5: Building elements from knots vector v = (0, 0, 0, 0.25, 0.50, 0.75, 1, 1, 1) 



Generally a open knots vector can be partitioned into rig elements expressed by 

ne = |E| = n - 2{k - 1) (94) 

4.7.2 Numerical Scheme for Time Integration 

Let the D-NURBS evolution equation: 

Mp + Dp + Kp = fp{p)- Ip (95) 

where /p = / J'^f{x, y, z)du and I{p) = f fiJ^Jdu, according to Appendix A and sections I4!4ll4.5l 
Let us consider the following numerical scheme: 

p(t+Ai) _ 2p(t) + p(t-At) 

P = 9 (96) 

(At? 



p{t+At) _ p(t-At) 

2Ai 



(97) 



If we substitute these expressions in equation (1951 ) we obtain: 



„(t+At) 9„(t) _|_ „(t-At) \ f ^(t+At) ^{t-At)\ 
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(98) 

where M, D, K and / are supposed to be computed at time t + At. We shall be careful about the term 

-')-^p( 

2At 



J I p(t+^t)p(t At) \ pQ^Q^^jj^g definition in expression (11151 ) and expression (|97] ) we can write: 



2At 



Therefore, we can rewrite equation (|99) as: 



2At 



du. 



By using the fact that j(*+At)p(t+At) — q simplify expression (llOOl i to: 
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Using the approximation: 
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we get: 
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du 



= (|/^ (J^)^*^^*^ jit+At)^it-At) _l^^jTYt+At) ^(i_At)p(i-At)^ 

So, by substituting expression (IIOII ) in (|98] ) it renders: 

^ j^ p(*+A.)-2pW+p(^-^*) ^| ^ ^ ^ p(t+A.)_p(^-A.) \ ^ ^^^^^^^^ ^ 



{Aty 



2At 



If we multiply both sides of expression ( |98] ) to x4 (At)^ and rearrange the terms we get: 



(101) 



(102) 



u 
(103) 



Um + 2AtD + 4 {Atf p(*+^*) = 4 {Atf fp + 8MpW - (3M - 2AtZ)) p(*-^*) - J ^J^c^^-^^^d 
This expression can be written as: 

^(t+Ai)p(t+Ai)^^(M-A<)^ (104) 



where: 



^(^t+At) ^ 4^ _^ 2/\tD + 4 (At)^ iT, (105) 



and, 



Therefore, once initial conditions p(0) = po and p(0) = vq are given, we can use the approximation: 

p(0) - p(0 - At) 



At 

to write: 



VQ. (107) 



p(-At) = p(0) - Atvo, (108) 

and, consequently, we can start the iterative scheme given by expression (I104I ). 

The complexity for computing the expression (I104l i depends on the algorithm for calculating the matrices M, 
D, K and the method used to solve the linear system. Considering that n is the number of control points, is 
number of elements, k is the polynomial order of NURBS basis and rig is the number of quadrature points, the 
algorithm implemented to compute the matrices M, D, K performs the following steps: 

1 . For e = 1 ... rig do 

(a) Compute the Jacobian matrix block for element "e" (complexity 0{ng * n * k)). 

(b) Compute mass matrix block for element "e" (complexity 0{ng * k"^)). 

(c) Compute damping matrix block for element "e" (complexity 0{ng * k"^)). 

(d) Compute stiffness matrix block for element "e" (complexity 0{ng * fc^)). 

Therefore, the asymptotic complexity of the whole algorithm is given by: 

O {ueUg {O {nk) + O {k^) + O + O [k^])) = O {neUgnk) . (109) 

We highlight that the computational cost of the D-NURBS evolution must also consider the numerical method 
for solving the linear system [T04] To compute 1 104 1 we have used conjugate gradient method whose complexity is 
O (n). Hence, we can conclude that the expression 1 104l has final computational complexity equal to O {neUgnk). 

5 Experimental Results 

We have developed an experimental environment based on the D-NURBS approach with constraints. In our setting 
we consider the case of an elastic wire with 10m length with negligible transverse section fixed at the ends. 

The NURBS curve geometry is instantiated using an open knot vector v = (0, 0, 0, 0, 0.25, 0.50, 0.75, 1,1,1,1) 
with basis functions of order = 4 (degree A; — 1 = 3). Therefore, following section |2l the spline space has di- 
mension n — fc + l = 10 — 4+1 = 7, which means that we have seven controls points. Each point of a NURBS 
curve is influenced by k control points. Therefore, to set geometric constraints that keep the wire fixed at the ends, 
we must let — 1 fixed control points at the ends of the curve. This can be cast in the linear constraint framework 
for D-NURBS developed in section [ 
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Table 1: Initial configuration of the control points and weights (generalized coordinates) for wire simulation: 

p(0) = {xo,yo,zo,wo;xi,yi,zi,wi; • • X6,yQ, Z6,wq). 

Iteration: 0[0.00%] 



5 



10 

Figure 6: Intial D-NURBS setup for simulation of the elastic wire fixed at the ends. 



Besides, we consider that the wire is subject to a gravitational field with value g = dSm/s"^ and define 
control points position and weights at t = according to tabled] Besides, we set p(0) = to complete the initial 
conditions for time integration. 



Figure |6] demonstrates the environment at time t = Os. Here physics parameters were defined as: a = 35, 
/? = 10, ^ = 30, 7 = 0. To perform spatial and time integration we define 10 points in Gauss quadrature and 
At = 0.008s, respectively. 



In our experiments we observed that the evolution of the weights Wi may cause unrealistic behaviors and 
instabiUty, as observed in Figure |7](a). As mentioned in [ Terzopoulos and Qin 1994] , the weights Wi may not have 
arbitrary finite real values. Negative values may vanish the denominator of the rational functions in expression |7] 
Besides, small weights values may lower the deformation energy [ Terzopoulos and Qin 1994) . Therefore, some 
constraint must be included in order to enforce some control in the weight vector evolution. 

In this work we implement this task by a very simple strategy: the generalized coordinate vector is updated 
by solving the expression l90l but the weight vector is always returned to its initial value; that means, Wi = I for 
i = 0, 1, • • •, 6, following Table [T] As observed in Figure |7](b), the wire evolution becomes (visually) acceptable 
in this case. 
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(a) Wire configuration at iteration t = 1, 10, 22 without constrain thie weights evolution. 



iteration: 1 [0.28%] 



iteration: 10 [2.78%] 



iteration: 22 [6.11%] 




(b) System configuration at iteration t = 1, 10, 22 when enforcing weights Wi = 1 after each iteration. 
Figure 7: D-NURBS behavior for unconstrained and constrained weight vector evolution. 



To study the dynamic evolution of the D-NURBS curve we choose a point in the center of the wire and 
followed its amplitude evolution in time. Figure [8a]shows its dynamic evolution without the presence of damping, 
while Figure [8b] illustrates the dynamic evolution with damping, where 7 = 5. As expected, the former reports 
a periodic evolution once there are not dissipative forces and the latter pictures an attenuation of the amplitude 
along the time due to the damping. 



To analyze the effects of the elasticity and stiffness of D-NURBS curve we increasing each parameter sep- 
arately. First, we leave j3 (see equation ISTI ) with the same value of the initial configuration and modify a. The 
Figures |9a] and |9b] show the results. 

Similarly, we modify /3 while a remains unchanged. The results are shown in figures ^and ^ We observe 
that the system is more sensitive respect to the parameter /3 than the parameter a. In fact, when increasing the 
parameter /3 from 11 to 15 we observe a drastic change in the amplitude evolution as highlighted when comparing 
Figures [9c] and [9dl On the other hand, when changing a from 70 to 105 (Figures [9a] and [9b] respectively) we did 
not observe a similar behavior. 



The expression (1109 l l shows that the number of control points has a fundamental role in the computational 
cost of the D-NURBS algorithm. Therefore, we perform a runtime analysis of D-NURBS evolution for different 
number of controls points. We set parameters to: k = 3, a = 35, /3 = 50, /i = 30, 7 = 1 and five points in Gauss 
quadrature. The host is a Intel Core I5-3210M at 2.5 Ghz, with 6 GB RAM running a Windows 7 (64bit). 

We take 360 iterations for each configuration and measure the corresponding CPU time. In order to compare 
the complexity given by expression (1109 l l and the CPU time for each simulation, we compute the following rates: 

rrj + l _ rpj 

P' = , (110) 
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(a) Amplitude evolution for D-NURBS withiout damping. 
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(b) Amplitude evolution of D-NURBS with damping. 
Figure 8: Amplitude evolution of elastic wire represented by D-NURBS. 
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(b) a = 105 and /3 = 10. 
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(d) a = 35 and ^ = 15. 



Figure 9: Sensitivity of D-NURBS amplitude respect to elasticity a and stiffness /3. 
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Figure 10: Rates for computational complexity an CPU time given by expressions (lllll ) and (II lOI ). respectively. 
The former is pictured by the blue plot while the latter by the red one. 

where is the CPU time for configuration j and is the asymptotic complexity for the same configuration; 
that means: 

= 360 * ninin^ . (112) 

The configuration (j = 1) has 20 control points and k = A and 360 D -NURBS iterations are performed. Next, 
for j = 2, we increase the number of control points by 5, keep k = A and perform 360 iterations of the algorithm 
again, and so on. The result is pictured on Figure [TO] where the dot blue curve shows the evolution of expression 
(lllll ) and the red line shows the evolution of expression (II 101 ). both for j = 1, 2, • • •, 28 (number of control points 
20, 25, . . ., 140). 



By observing Figure [TO] we note that as we increase the control points number we get — )• . This is the 
expected behavior for asymptotic function, i.e., for large n the similarity between real and predicted time becomes 
more evident. 

6 Conclusions and Future Works 

We present a review of D-NURBS approach. We emphasize the formulation based on the Lagrangian mechanics 
followed by detailed development of the governing equations. We used a numerical method based on isogeo- 
metric analysis for the spatial integration used to compute the Jacobian, mass, damping and stiffness matrix. For 
validation we performed experiments with D-NURBS curve and discuss the influence of parameters, effects of 
NURBS weights and computational cost. For further works we plan to evaluate the D-NURBS for 2D and 3D 
systems. 
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A Appendix 

Let us consider the expression: 




(113) 



By the product rule we have: 




(114) 



Now, let us define the expressions / and / as: 
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Therefore, we can rewrite expression (1113^ : 
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Now, we prove that: 
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If we name: 
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then, by applying the product rale we get: 
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(120) 



Let jj be the coUum i of the Jacobian J. Once J = J (p) , where p = p (t) , then the Chain-Rule allows to 
write: 



d ^.^ dJ d ^ ^ dJ ■ 

dt dpi dt dpi 



Expression (11201 ) can be rewritten as: 



I f dJ ■ 



T dJ ■ 



dp, 



P, 



So, by substituting equation (11211 ) in expression in (11221 ) we obtain: 



R=l{Mfjp + l{jpfi3 



I ] 1 
T 



^ ^pY ^{J''J)P^ /or i = 0,...,4(n + l), 



Therefore, from the expressions (11 19l l and (I124t we get that 

{hfjp=^^[p) 

which is equivalent to expression (II 181 ). 

Therefore, by substitution this result in equation (I117l i we obtain: 

Y{p,p) = Ip, 

where / is computed by expression (I115I ). 
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B Appendix 

In order to prove that: 



we must observe that: 
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However, due to the definition of jacobian matrix J we must have ^ = ji- Therefore: 

dJ 



dp., 



p = 0. 



On the other hand: 
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But, from the definition of K matrix in expression (1561 ): 
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Once the vector p does not depend on the parameter u we can write the first term inside the integral as: 

tT 



due to equation (I129I ). Therefore, expression (I127l i has been proved. 



